  function [v,s,perc,sim]=see_igone( mu , sd , lb , ub , ndr );
% function [v,s,perc,sim]=see_igone( mu , sd , lb , ub , nd );
% See IG1 with mean = MU 
% and for the time being infinite variance 
% lb  (optional) lower bound of the graph (default = 0.01)  
% ub  (optional) upper bound of the graph (default = 5 )
% ndr (optional) number of draws ( default = 20000 ) 
% V=2 can change that from loc_gone.m 
% Alejandro Justiniano
% Last Modified December 7 2006 
% ========================================================
pvec=[0.005 0.01 0.05 0.1 0.5 0.9 0.95 0.99 0.995]; 
if nargin < 5 
    ndr=100000; 
    if nargin < 4 
        ub = 8 ;
        if nargin < 3 
            lb = 0.0001; 
        end
    end
end
if sd > 40 
    sd=inf 
end 
[svec,vvec]=inverse_gamma_specification( mu , sd ,1 ); 
gr=linspace(lb,ub);
jj=1; 
for jj=1:length(vvec); 
    v=vvec(jj); 
    s=svec(jj); 
    pd=zeros(100,1); 
    ii=1; 
    for ii=1:100; 
        pd(ii)=pdf_igone(gr(ii),v,s) ; 
    end 
    mod=( (v/(v+1))^0.5 )*sqrt(s/v); 
    
    into.title=['IG1 Mean ',num2str(mu),' SD: ',num2str(sd), ' ; v ', num2str(v) , ' s= ',num2str(s),' mode : ',num2str(mod) ]; 
    subplot(2,1,1);
    plot(gr,pd); 
    title([' PDF ',into.title ] );
    dr=zeros(ndr,1); 
    sim=1./gamrnd(v/2,2/s,ndr,1); 
    sim=sqrt(sim);
    display( ['Mean Simulated ' num2str( mean(sim) ) '; Std simulated ' num2str( std(sim) ) ] ) ;
    pvec=floor(ndr*pvec); 
    sims=sort(sim); 
    perc=[(pvec/ndr)' sims(pvec)]; 
    nper=length(pvec); 
    cper=emptycell(nper,1); 
    
    std_sim = std( sim ) 
    
    
    subplot(2,1,2); 
    %binedge=(lb:((ub-lb)/50):ub); 
    %[N,BIN]=histc(sim,binedge);
    %N=N/( length(sim) );
    %bar(binedge,N,'histc');
    temp1=min(find( sims > lb ));
    temp2=max(find( sims < ub )); 
    
   
    hist( sims(temp1:temp2) , 100 ); 
    ii=1;
    for ii=1:nper;
        temp1=sprintf('%3.3f',perc(ii,1) );
        temp2=sprintf('%3.3f',perc(ii,2) );
        cper(ii)={ [temp1,'  ',temp2] };
    end;
    legend( cper{1}, cper{2}, cper{3},  cper{4}, cper{5}, cper{6}, cper{7}, cper{8},  cper{9} );
    %histplot(sim , , into );
    printcell(cper); 
end;
    
